%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Program: tables.m
% By: Stephie Fried, Kevin Novan, and William Peterman
% Date: Summer 2022
% Purpose: Produces tables and figures for paper
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
close all;
clear all;

bpath = 'C:\Users\l1sxf17\Dropbox (FRB SF)\Research\Will\equity_efficiency_project\Replication';
figpath ='C:\Users\l1sxf17\Dropbox (FRB SF)\Research\Will\equity_efficiency_project\Replication\Figures_tables';
tablepath ='C:\Users\l1sxf17\Dropbox (FRB SF)\Research\Will\equity_efficiency_project\Replication\Figures_tables';

pic = 1; 
cd(bpath)
%% Table 1

%Load processed results
load ('Programs/Matlab/aggregates')

cd(tablepath);
FID = fopen('table1.tex', 'w');
fprintf(FID, '\\begin{table}[H] \n');
fprintf(FID, '\\center \n');
fprintf(FID, '\\singlespace \n');
fprintf(FID, '\\caption{Distribution and Welfare Effects} \\label{cev}\n');
fprintf(FID, '\\vspace{.05in} \n');
fprintf(FID, '\\begin{tabular}{l c c c c c}\\hline \n');
fprintf(FID, '\\hline \n');
fprintf(FID, ' & \\multirow{2}{*}{CEV}  & \\multicolumn{4}{c}{Percent change in}\\\\ [-.5ex] \n');
fprintf(FID, '\\cline{3-6} \n');
fprintf(FID, ' & & Welfare Gini & Output & Capital & Energy \\\\ [.5ex] \n');
fprintf(FID, '\\hline \n');
fprintf(FID, 'Optimal &  %8.2f &  %8.1f  &  %8.1f & %8.1f & %8.1f\\\\ ', ...
    policies(cev_ind, 1), (policies(giniw_ind,1)./giniw_base -1)*100, (output(1)/Y_base-1)*100,...
     (policies(K_ind,1)./K_base -1)*100,  (policies(E_ind,1)./E_base -1)*100);
 fprintf(FID, 'Capital tax &  %8.2f &  %8.1f  &  %8.1f & %8.1f & %8.1f\\\\ ', ...
    policies(cev_ind, 2), (policies(giniw_ind,2)./giniw_base -1)*100, (output(2)/Y_base-1)*100,...
     (policies(K_ind,2)./K_base -1)*100,  (policies(E_ind,2)./E_base -1)*100);
 fprintf(FID, 'Labor tax &  %8.2f &  %8.1f  &  %8.1f & %8.1f & %8.1f\\\\ ', ...
    policies(cev_ind, 3), (policies(giniw_ind,3)./giniw_base -1)*100, (output(3)/Y_base-1)*100,...
     (policies(K_ind,3)./K_base -1)*100,  (policies(E_ind,3)./E_base -1)*100);
 fprintf(FID, 'Lump sum &  %8.2f &  %8.1f  &  %8.1f & %8.1f & %8.1f\\\\ ', ...
    policies(cev_ind, 4), (policies(giniw_ind,4)./giniw_base -1)*100, (output(4)/Y_base-1)*100,...
     (policies(K_ind,4)./K_base -1)*100,  (policies(E_ind,4)./E_base -1)*100);
 fprintf(FID, 'Labor progressivity &  %8.2f &  %8.1f  &  %8.1f & %8.1f & %8.1f\\\\ ', ...
    policies(cev_ind, 5), (policies(giniw_ind,5)./giniw_base -1)*100, (output(5)/Y_base-1)*100,...
     (policies(K_ind,5)./K_base -1)*100,  (policies(E_ind,5)./E_base -1)*100);
 fprintf(FID, 'Income dependent &  %8.2f &  %8.1f  &  %8.1f & %8.1f & %8.1f\\\\ ', ...
     policies(cev_ind, 6), (policies(giniw_ind,6)./giniw_base -1)*100, (output(6)/Y_base-1)*100,...
     (policies(K_ind,6)./K_base -1)*100,  (policies(E_ind,6)./E_base -1)*100);
fprintf(FID, '\\hline \n');
fprintf(FID, '\\end{tabular}\n');
fprintf(FID, '\\end{table} \n');
fclose(FID);
cd(bpath);


%% Numbers for Table 3: Sensitivity

cd(tablepath);
FID = fopen('table3.tex', 'w');
fprintf(FID, '\\begin{table}[H] \n');
fprintf(FID, '\\center \n');
fprintf(FID, '\\singlespace \n');
fprintf(FID, '\\caption{Sensitivity} \\label{tab:sensitivity}\n');
fprintf(FID, '\\begin{tabular}{l c c}\\hline \n');
fprintf(FID, '\\hline \n');
fprintf(FID, ' & Fraction of revenue used  & Percent change in\\\\ \n');
fprintf(FID, ' &to reduce the capital tax & the welfare Gini\\\\ \n');
fprintf(FID, '\\hline \n');
fprintf(FID, ' \\emph{Subsistence energy: $\\bar{e}$} \\\\ \n');
fprintf(FID, '\\hline \n');
fprintf(FID, '$\\bar{e}=0$ &  %8.2f &  %8.2f\\\\ ',...
     (lambdak_base*r_base*K_base - EBARZERO(lambdak_ind)*EBARZERO(r_ind)*EBARZERO(K_ind))/(taue*EBARZERO(E_ind)),... 
      round((EBARZERO(giniw_ind)/giniw_base-1)*100, 2));
fprintf(FID, '$\\bar{e}=0.0013$ &  %8.2f &  %8.2f\\\\ ',...
     (lambdak_base*r_base*K_base - OPT(lambdak_ind)*OPT(r_ind)*OPT(K_ind))/(taue*OPT(E_ind)),... 
     round((OPT(giniw_ind)/giniw_base-1)*100, 2));  
 fprintf(FID, '$\\bar{e}=0.0026$ &  %8.2f &  %8.2f\\\\ ',...
     (lambdak_base*r_base*K_base - EBAR2X(lambdak_ind)*EBAR2X(r_ind)*EBAR2X(K_ind))/(taue*EBAR2X(E_ind)),... 
     round((EBAR2X(giniw_ind)/giniw_base-1)*100, 2));  
 fprintf(FID, '\\hline \n');
fprintf(FID, ' \\emph{Carbon tax: $\\tau^c$} \\\\ \n');
fprintf(FID, '\\hline \n');
fprintf(FID, '\\$30/ton CO$_2$  &  %8.2f &  %8.2f\\\\ ',...
     (lambdak_base*r_base*K_base - SMALLTAX(lambdak_ind)*SMALLTAX(r_ind)*SMALLTAX(K_ind))/(0.75*taue*SMALLTAX(E_ind)),... 
     round((SMALLTAX(giniw_ind)/giniw_base-1)*100, 2));  
 fprintf(FID, '\\$40/ton CO$_2$ &  %8.2f &  %8.2f\\\\ ',...
     (lambdak_base*r_base*K_base - OPT(lambdak_ind)*OPT(r_ind)*OPT(K_ind))/(taue*OPT(E_ind)),... 
     round((OPT(giniw_ind)/giniw_base-1)*100, 2));  
fprintf(FID, '\\$50/ton CO$_2$  &  %8.2f &  %8.2f\\\\ ',...
     (lambdak_base*r_base*K_base - BIGTAX(lambdak_ind)*BIGTAX(r_ind)*BIGTAX(K_ind))/(1.25*taue*BIGTAX(E_ind)),... 
     round((BIGTAX(giniw_ind)/giniw_base-1)*100, 2)); 
fprintf(FID, '\\hline \n');
fprintf(FID, '\\end{tabular}\n');
fprintf(FID, '\\end{table} \n');
fclose(FID);
cd(bpath);

%% Table 4: Fraction Better off

%Load fraction better off (produced in individual_cevs2.m)
load('Programs/Matlab/frac_better_off');
%Order for frac_better_off: {'optimal', 'simple_capital', 'simple_labor', 'simple_lump', 'optimal_prog', 'MEID'};

cd(tablepath);
FID = fopen('table4.tex', 'w');
fprintf(FID, '\\begin{table}[H] \n');
fprintf(FID, '\\center \n');
fprintf(FID, '\\singlespace \n');
fprintf(FID, '\\caption{Share of Agents Better Off} \\label{cev}\n');
fprintf(FID, '\\vspace{.05in} \n');
fprintf(FID, '\\begin{tabular}{l c c c c c c }\\hline \n');
fprintf(FID, '\\hline \n');
fprintf(FID, ' & Capital & Labor   &  Lump  & Labor   &   Income & \\multirow{2}{*}{Optimal} \\\\ [-.5ex] \n');
fprintf(FID, ' & tax & tax & sum &   progressivity  & dependent \\\\ [.5ex] \n');
fprintf(FID, '\\hline \n');
fprintf(FID, '\\%% better off  &  %8.1f &  %8.1f  &  %8.1f & %8.1f & %8.1f  & %8.1f\\\\  ',...
   frac_better_off(2:6)*100,  frac_better_off(1)*100);
fprintf(FID, '\\%% prefer optimal  &  %8.1f &  %8.1f  &  %8.1f & %8.1f & %8.1f  & -\\\\  ',...
   frac_prefer(2:6)*100);
fprintf(FID, '\\hline \n');
fprintf(FID, '\\end{tabular}\n');
fprintf(FID, '\\end{table} \n');
fclose(FID);
cd(bpath);

%% Figure 1

%Size of the progressive rebate
inc_ratio = [0.01:0.0005:3];
base_tax  = max(1 - lambda0_base*inc_ratio.^(-lambda1_base), 0);
optimal_tax =  max(1 - OPT(lambda0_ind)*inc_ratio.^(-OPT(lambda1_ind)), 0);
prog_rebate = max((base_tax - optimal_tax)./base_tax*100, 0);

%Kink point
ind1 = find(prog_rebate ==0, 1, 'first');
kink_point = inc_ratio(ind1);

avg_rate = 1 - lambda0_base*inc_ratio.^(-lambda1_base);
avg_rate_opt = 1 - OPT(lambda0_ind)*inc_ratio.^(-OPT(lambda1_ind));
avg_rate_opt(avg_rate_opt <0)=0;
avg_rate_opt(avg_rate_opt>avg_rate) = avg_rate(avg_rate_opt>avg_rate);

figure(1)
hax=axes; 
hold on
plot(inc_ratio(1:50:end), avg_rate(1:50:end), '-x', 'LineWidth', 2)
plot(inc_ratio(1:1:end), avg_rate_opt(1:1:end),  'LineWidth', 2)
xlabel('Labor Income Relative to Mean Labor Income')
xlim([0, 1.5])
ylabel({'Average labor tax rate'})
legend('Baseline', 'Optimal', 'Location', 'NW')
legend boxoff
set(gca, 'Fontsize', 12)
box off
if(pic ==1)
    saveas(gcf,[figpath ,filesep, 'avg_tax_rate'],'epsc');
end

%% Figure 2
load ('Programs/Matlab/mu')

figure(2)
hold on
plot(cumpop*100, cs_prog,'--','linewidth', 2, 'Color', [31,120,180]/255)
plot(cumpop*100, cs_idme, 'linewidth', 2, 'Color', [166,206,227]/255)
plot(cumpop*100, cs_lump, '-.', 'linewidth', 2, 'Color', [178,223,138]/255)
xlim([0, 100])
ylim([0, 1])
set(gca, 'Fontsize', 12)
ylabel('Cumulative Share of Revenue')
xlabel('Marginal Utility Percentile (ranked highest to lowest)')
legend('Labor-Progressivity Rebate', 'Income-Dependent Rebate',  'Lump-Sum Rebate', 'Location', 'SE')
legend boxoff
box off
if(pic ==1)
    saveas(gcf,[figpath ,filesep,'mu_figure'],'epsc');
end


%% Figure 3:  Ex-post distribution 
load ('Programs/Matlab/results_cevi'); 
mat = [welfare_base, cevimat*100, weldiff]; 
mat2 = sortrows(mat, 1);

welfare_base_sorted = mat2(:, 1); 
opt_sorted = mat2(:,2); 
cap_sorted = mat2(:, 3); 
labor_sorted  =mat2(:,4); 
lump_sorted = mat2(:, 5); 
prog_sorted = mat2(:, 6); 
meid_sorted = mat2(:, 7); 

np =100; %size of the bin
nb = ceil(10000/np);
opt_ile = zeros(nb, 1); 
labor_ile = zeros(nb, 1); 
lump_ile = zeros(nb, 1); 
cap_ile = zeros(nb, 1);
prog_ile = zeros(nb, 1); 
meid_ile = zeros(nb, 1); 

for i = 1:nb
    ind1 =(i-1)*np +1; 
    ind2 = min(i*np, 10000); 
    opt_ile(i) = mean(opt_sorted(ind1:ind2)); 
    cap_ile(i) = mean(cap_sorted(ind1:ind2)); 
    labor_ile(i) = mean(labor_sorted(ind1:ind2)); 
    lump_ile(i) = mean(lump_sorted(ind1:ind2));
    prog_ile(i) = mean(prog_sorted(ind1:ind2)); 
    meid_ile(i) = mean(meid_sorted(ind1:ind2)); 
end

figure(3)
hold on
plot(1:nb, opt_ile, '-k', 'LineWidth', 1.5)
plot(1:nb, cap_ile,  ':', 'LineWidth', 1.5, 'Color',[0.4660, 0.6740, 0.1880])
plot(1:nb, labor_ile, '-.', 'LineWidth', 1.5, 'Color', 	[0.3010, 0.7450, 0.9330])
xlabel('Lifetime Welfare Percentile')
ylabel('Percent')
title('Efficiency Instruments')
legend('Optimal', 'Capital tax', 'Labor tax')
legend boxoff
set(gca, 'Fontsize', 12)
if(pic ==1)
    saveas(gcf,[figpath ,filesep, 'figure2a'],'epsc');
end

figure(4)
hold on
plot(1:nb, opt_ile, '-k', 'LineWidth', 1.5)
plot(1:nb, lump_ile, '--','LineWidth', 1.5)
plot(1:3:nb, prog_ile(1:3:end), '-x','LineWidth', 1)
plot(1:3:nb, meid_ile(1:3:end), '-*','LineWidth', 1)
xlabel('Lifetime Welfare Percentile')
title('Equity Instruments')
ylabel('Percent')
legend('Optimal', 'Lump sum', 'Labor progressivity' ,'Income dependent')
legend boxoff
set(gca, 'Fontsize', 12)
if(pic ==1)
    saveas(gcf,[figpath ,filesep, 'figure2b'],'epsc');
end

%% Figure 4: Transition 

%%% Lifetime C
M = xlsread("Fortran_output/transition/summary of results", 'constant C');
age = M(:,2); 
opt_trans_cev = M(:, 3); 
combo_trans_cev = M(:,5); 

age1 = [-80:1:80]; 
figure(5)
hold on
plot(age1, [opt_trans_cev(1)*ones(size(age1(1:41)))'; opt_trans_cev], 'linewidth', 2, 'Color', [166,206,227]/255)
plot(age1,  [opt_trans_cev(1)*ones(size(age1(1:41)))'; combo_trans_cev], '--','linewidth', 2, 'Color',  [31,120,180]/255)
plot(age1, policies(cev_ind, 1)*ones(size(age1)),'-.','linewidth', 1, 'Color', [0,0,0]/255)
legend('Optimal: transition', 'Phase-in: transition', 'Steady state', 'Location', 'NW')
legend boxoff
xlabel('Model age at time of adoption')
ylabel('Lifetime CEV (percent)')
set(gca, 'Fontsize', 12)
ylim([-0.3, 0.31])
if(pic ==1)
    saveas(gcf,[figpath ,filesep, 'trans'],'epsc');
end

